function [Gamma_g,Su_g]=BVARgibbs_draws(Par,nGibbs)

CSuInv=Par.CSuInv;
CZZinv=Par.CZZinv;
Gamma=Par.Gamma;
T=Par.T;
Tpr=Par.Tpr;
k=Par.k;
N=Par.N;
p=Par.p;
jg=1;

while jg <= nGibbs
    
    U    = randn(T+Tpr-k+2,N)*CSuInv;                % Inverted Wishart and Variance
    Su_temp  = inv(U'*U);
        
    temp  = randn(size(Gamma));                       % Compute the betas
    Gamma_temp = Gamma + CZZinv'*temp*chol(Su_temp);  % Instead of using KRONECHER PRODUCT!!!!!!
%     
%     AA=zeros(N*(p+1));
%     AA(1:N,1:N*p)=Gamma_temp(1:end-1,:)';
%     AA(N+1:N*p,1:N*(p-1)) = eye(N*(p-1));
%     AA(end-N+1:end,end-N+1:end)=eye(N);
%     AA(1:N,end-N+1:end)=eye(N);
    Gamma_g(:,:,jg) = Gamma_temp;
    Su_g(:,:,jg)=Su_temp;
    jg = jg+1;
    
end;

%Su_g = Su_g(:,:,1:end);
%Gamma_g = Gamma_g(:,:,1:end);